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We present a model supported by simulation to explain the effect of temperature on the conduction 
threshold in disordered systems. Arrays with randomly distributed local thresholds for conduction 
occur in systems ranging from superconductors to metal nanocrystal arrays. Thermal fluctuations 
provide the energy to overcome some of the local thresholds, effectively erasing them as far as the 
global conduction threshold for the array is concerned. We augment this thermal energy reasoning 
with percolation theory to predict the temperature at which the global threshold reaches zero. We 
also study the effect of capacitive nearest-neighbor interactions on the effective charging energy. 
Finally, we present results from Monte Carlo simulations that find the lowest-cost path across an 
array as a function of temperature. The main result of the paper is the linear decrease of conduction 
threshold with increasing temperature: Vt{T) = Vt(0)(l - A.SkBTP{Q)/pc), where 1/P(0) is an 
effective charging energy that depends on the particle radius and interparticle distance, and Pc is 
the percolation threshold of the underlying lattice. The predictions of this theory compare well to 
experiments in one- and two-dimensional systems. 
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I. INTRODUCTION 

In many physical systems, local barriers prevent the 
onset of steady-state motion or conduction unless a cer- 
tain minimum threshold for an externally applied driving 
force or bias is exceeded. Often, the strength of those 
barriers varies throughout the system and only their sta- 
tistical distribution is known. A key issue then concerns 
how the global threshold for onset of motion is related to 
the distribution of local threshold values. Examples in- 
clude the onset of resistance due to depinning of fluxline 
motion in type-II superconductors, the onset of mechan- 
ical motion in coupled frictional systems such as sand 
piles, and the onset of current flow through networks of 
tunnel junctions in the Coulomb blockade regime. In all 
of these cases, defects in the host material or the under- 
lying substrate produce local traps or barriers of varying 
strength. 

Under an applied driving force, fluxlines, mobile parti- 
cles or charge carriers from an external reservoir can pen- 
etrate the disordered energy landscape, becoming stuck 
at the traps or piling up in front of barriers. With in- 
creased drive, particles can surmount some of the barri- 
ers and penetrate further. However, a steady-state flow 
is only established once there is at least one continuous 
path connecting one side of the system with the other. 
The onset of steady-state transport then corresponds to 
finding the lowest-energy system-spanning path. This 
optimization problem was addressed in 1993 in a seminal 
paper by Middleton and Wingren (MW)A 

Using analytical arguments as well as computer simu- 
lations, MW found that, for the limit of negligible ther- 
mal energies, the onset of system-spanning motion corre- 
sponds to a second order phase transition as a function of 
applied bias. The global threshold value scales with dis- 
tance across the system, but is independent of the details 



of the barrier size distribution. Beyond threshold, more 
paths open up and the overall transport current increases. 
As a result, the steady-state transport current displays 
power law scaling as a function of excess bias. These 
predictions have subsequently been used extensively in 
the interpretation of single electron tunneling data from 
networks of lithographically defined junction arraysSi^ as 
well as from self- assembled nanoparticle systems In 
addition, recent experiments^ and simulations^ have ex- 
plored how the power law scaling is affected by structural 
disorder in the arrays. The regime of large structural dis- 
order and significant voids in the array was investigated 
numerically using a percolation model >2i 

What happens at finite temperature? Intuitively, one 
might expect temperature to produce a smearing of the 
local thresholds and thus a quick demise of the power 
law scaling for T > 0. Indeed, a number of experiments 
have found that the nonlinear current-voltage charac- 
teristics observed at the lowest temperatures give way 
to nearly linear, Ohmic behavior once T is raised to a 
few dozen Kelvin.^'''-'^''- More recently, however, several 
experiments showed that the scaling behavior survives 
with a well-defined, albeit now temperature-dependent, 
global threshold. In a previous Letter, we demonstrated 
for a two-dimensional metal nanocrystal array that a) 
the threshold is only weakly temperature dependent, de- 
creasing linearly with increasing T, and b) the scaling 
exponent remains unaffected by temperature. Conse- 
quently, the shape of the nonlinear response as a function 
of applied drive remains constant and is merely shifted 
to lower drive values as T increased 

Similar behavior was also observed in small 2D metal 
nanoparticle networks by Ancona et al~ and Cordan et 
al}^ and in ID chains of carbon particles by Bezryadin. 
et alm^ Most recently, it was corroborated by simula- 
tions of (semi-classical) particles in 2D arrays of pinning 



2 



sites with random strengthssi^ This weak temperature 
dependence of the nonhnear response also has important 
practical consequences as it implies that arrays are much 
more robust and forgiving as compared to systems with 
a single threshold that might be significantly affected by 
its local environment. 

However, the theoretical approach developed by MW 
considers only the zero-temperature limit where the lo- 
cal energy levels are sharply delineated and barriers be- 
tween adjacent sites are well-defined. In Ref. we in- 
troduced the main results from a new model that extends 
the MW approach to finite temperatures. Here, we de- 
velop this model in more detail, providing both analytical 
results and data from computer simulations. For con- 
creteness, we focus on single electron tunneling through 
metal nanoparticle arrays. However, we expect the main 
results to carry over to a much wider class of systems 
with distributed thresholds due to quenched disorder. 

Our model goes beyond previous work in two impor- 
tant aspects. First, we introduce a method that allows 
us to treat finite temperatures. This method is based on 
estimating when small barriers, washed out by tempera- 
ture, have percolated across the system, and it establishes 
an upper limit on the global threshold as a function of 
temperature. A key finding is that random quenched dis- 
order leads to universal behavior that is independent of 
the details of the barrier height distribution. Second, we 
include nearest neighbor capacitive coupling. This leads 
us to a new definition of the relevant effective charging 
energy for system crossing, in terms of the most proba- 
ble value in the distribution of energy costs. As shown in 
Ref. 0, the model captures the experimentally observed 
temperature dependence of the drive-response character- 
istics and predicts the collapse of the global threshold as 
a function of temperature on a universal curve that is 
independent of local junction parameters. 

The paper is organized as follows. In Section II we out- 
line the basic ingredients of our model, mainly focusing 
on the limit of negligible interparticle coupling. Section 
III then calculates the shape of the probability distribu- 
tion of energy costs for the general case of finite nearest 
neighbor coupling. In Section IV we present simulation 
results for various network geometries. We also discuss 
the validity of the percolation model and show numerical 
results for the decrease of the threshold with tempera- 
ture. Section V describes how the current-voltage char- 
acteristics behave at temperatures above the point where 
the voltage threshold reaches zero. Section VI contains 
a discussion of the model and comparisons with recent 
experimental data and as well as with numerical results 
from related systems. 



II. THE BASIC MODEL IN THE ABSENCE OF 
INTERPARTICLE CAPACITIVE COUPLING 

We consider one- or two-dimensional arrays of spheri- 
cal metal nanoparticles ("sites"), placed between two in- 



plane metal electrodes. We ignore any particle-internal 
level spacing due to quantum size effects and treat each 
site as possessing a continuous spectrum of available 
states up to some local chemical potential. This is a rea- 
sonable approximation for metal particles with diameters 
larger than a couple nanometers at temperatures above 
liquid helium. For such particles, the largest energy be- 
sides thermal energy is the electrostatic energy associated 
with the transfer of additional, single electrons. 

We consider interparticle spacings small enough to al- 
low for such transfer by electron tunneling. We make the 
usual assumptions of the "orthodox theory" of single elec- 
tron tunneling (see, e.g., Likharev in Ref. Il6|) . namely 
that the tunnel time is negligible in comparison with all 
other time scales, the tunnel resistance R >> Rq = /i/e^, 
where Rq is the quantum of resistance, co-tunnel events 
due to coherent quantum processes can be ignored, and 
the local tunnel rate from site to site depends only on 
the change in electrostatic free energy of the system, AE, 
that would result from a tunnel event. At low tempera- 
tures, a positive AE implies a suppression of tunneling 
(Coulomb blockade), and current flows only after an ex- 
ternal bias has been applied that compensates for this 
energy cost. If tunneling occurs from a site at higher en- 
ergy to one at lower energy {AE < 0), we assume that 
the energy difference is lost due to scattering processes 
in the destination particle (inelastic tunneling). 

Throughout the paper, we consider the limit of negli- 
gible structural disorder of the arrays, i.e., all sites are 
identical in terms of both their tunnel coupling and ca- 
pacitive coupling to neighbors, as well as in their self- 
capacitance. Disorder enters in form of a random distri- 
bution of the local chemical potentials at every site due 
to quenched offset charges. This quenched charge disor- 
der models charge fluctuations due to impurities in the 
substrate which in turn polarize the nanoparticles. 

A corresponding experimental system can be realized 
as shown in Ref. by self-assembling, onto an insu- 
lating substrate, ligand-coated nanoparticles from solu- 
tion. The ligands prevent nanoparticle sintering and 
well-ordered arrays are formed through a balance be- 
tween attractive van der Waals forces and repulsive steric 
hindrance between ligands from neighboring particles. 
For dodecancthiol ligands and particle diameters in the 
range 4.5nm to 7nm, a size dispersion of less than 5% 
can be achieved, resulting in 2D arrays with excellent 
long-range order of the particle packing. Electronic mea- 
surements, both on nanoparticle arrays but also on self- 
assembled monolayers of molecules by themselves, have 
shown that alkanethiol ligands act as mechanical spacers 
and do not otherwise affect the transport properties li^ii^ 
Consequently, they set the width of the tunnel barrier 
between neighboring nanoparticles but do not introduce 
states inside the barrier. 

The quenched charge disorder is not a perturbative ef- 
fect: in principle, the chemical potential of a nanoscale 
particle can be shifted by a nearby trapped charge as 
much as it would be by an added mobile electron. There- 
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fore, electrons in an array propagate through a network 
of junctions with randomly varying threshold voltages. 
Note that the mobile charges are quantized (electrons) 
and thus move the local chemical potential by the same 
amount, A/i, every time a single charge enters or leaves 
a site. On the other hand, the quenched charges model a 
polarization effect and thus can move the local chemical 
potentials continuously, just like a local gate electrode 
could. The overall energy cost, AE, associated with a 
tunnel event therefore has to take into account the ef- 
fect of both discrete mobile charges and of a continuous 
random distribution of quenched charges. 

One might expect conduction through large arrays to 
depend on the details of the local quenched, or back- 
ground, charge distribution. However, zero-temperature 
arguments by MW indicate that this is not the casCji as 
least in the limit of negligible capacitive coupling between 
sites. Instead, the overall array current-voltage charac- 
teristics (IVs) appear to be robust to background charge 
disorder and exhibit a non-zero effective voltage thresh- 
old, Vj, that scales linearly with array size (i.e., distance 
between electrodes). 

To see this, consider first a ID array at T = with 
a given distribution of quenched polarization charge val- 
ues. Because mobile electrons can compensate for lo- 
cal polarizations in integer multiples of e, the electronic 
charge, only disorder in the range [— e/2, +e/2] needs to 
be considered. Starting from an initial state of zero ap- 
plied overall bias, mobile charges can penetrate, say from 
the left, single-file into the disordered potential landscape 
until they first encounter a local up-step in electrostatic 
potential, AV > 0. At this point, the Coulomb blockade 
prevents further advance. 

To the left of the up-step, each site now has one addi- 
tional charge on it and all potentials have been raised uni- 
formly by e/Co where Co is the self-capacitance of each 
site. In order to move the charge front further toward 
the right electrode, the bias applied to the left electrode 
has to be raised. Each time an up-step is encountered 
anywhere in the array, a bias increment of e/Cg at the 
left electrode will suffice to advance the front. Thus, the 
minimum bias required in order for mobile charges to 
make it all the way across the array will be given by the 
number of up-steps times e/Co (recall that down-steps in 
local potential do not matter as tunneling is assumed to 
be inelastic). In other words, the T = 0, global threshold 
for conduction for an array of N sites is given by Ref. Q 
as 

Vt{0) ^ aNe/Co. (1) 

If we now assume a flat, random distribution of quenched 
charges, on average half of the steps between neighboring 
sites will be up-steps. Therefore, for ID arrays a = 1/2. 

Note that this argument of MW depends only on the 
number of up-steps, but not on their magnitude |AV^|! 
Thus, details of the distribution of step sizes are irrele- 
vant at T = 0. This also holds for 2D systems, except 
that now the mobile charges can, to some extent, avoid 



up-steps. Consequently, there will be some roughness 
in the front of charges advancing across the array below 
threshold. Equation ^ still holds, with N now the dis- 
tance across the gap between the electrodes. Na is the 
number of up-steps in the path across the array with the 
least number of up-steps ( "optimal path" ) . The value of 
a in 2D will be smaller than in ID and depend on the ar- 
ray topology. Unfortunately, analytical arguments that 
would predict a for 2D systems are not known and one 
has to resort to computer simulations. Specifically, for 
a close-packed triangular arrangement of spheres we find 
a = 0.226 (see Section IV). 

In order to model the effect of finite temperature on 
the global threshold for conduction, we start by consid- 
ering thermal fluctuations at the local, single junction 
level. Let AE denote the change in the electrostatic po- 
tential energy of the system when a single electron moves 
from one site to another. If | Ai?| >> fc^T, the nonlinear. 
Coulomb blockage dominated current-voltage character- 
istic will survive: current will be suppressed below the 
local voltage threshold but will rise approximately lin- 
early above it.^*^ On the other hand, for |A_E| << fc^T, 
the Coulomb blockade vanishes and the junction conduc- 
tance will exhibit linear, Ohmic behavior down to the 
lowest bias voltage. 

As a first approximation, we now coarse-grain the sys- 
tem into two categories of tunnel junctions. Junctions 
between sites with energy differences |Ai?| > bksT will 
be treated as if T = 0, implying a fully nonlinear re- 
sponse and, below threshold, the absence of zero-bias 
conductance. Junctions between sites with energy dif- 
ferences \AE\ < hksT will be treated as if AE — and 
all Coulomb blockade effects were removed, implying a 
linear response like Ohmic conductors. The parameter h 
measures the extent of thermal broadening and depends 
on details of the electronic level distribution. If energy 
levels are within bksT, then electrons from thermally 
excited states above the Fermi level on site i can tun- 
nel directly into available states below the Fermi level 
on neighboring site j. This means that up-steps within 
hksT are effectively removed. 

To determine 6, we consider in each nanoparticle the 
width of the tail of unoccupied states below and of oc- 
cupied states above the Fermi level. Each tail has an 
approximate width of ksT so that \AE\ is reduced by 
roughly iksT and thus & w 2. To make this argument 
more quantitative, we consider the mean energy of states 
above the Fermi energy /i in particle i, 

^ j^EimmdE 

^^'"sh}^- j:)^E)f{E)dE 

where D{E) is the density of states and f{E) is the 
Fermi-Dirac function. Evaluating the integral as a series 
and determining the coefficients numerically, we obtain 
{Ehigh)i ~ fJ-i + l.iksT. By symmetry, the mean energy 
of the low-energy unoccupied tail in particle j will be 
{Eiow)j ~ — 1.2kBT. Tunneling from the high-energy 
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tail of particle i to the low-energy tail of particle j thus 
will cost a mean energy AE — {fij — fii) — 2AkBT. This 
leads to 5 = 2.4. 

As temperature is raised, more and more junctions will 
satisfy |Ai?| < bksT and lose their nonlinear behavior. 
We define p{T) as the fraction of junctions that has been 
effectively linearized. Since both up- and down-steps will 
be affected equally by thermal smearing, p{T) can be 
found from 

pbkBT 

p(T) = 2 / P{AE)dAE (2) 
Jo 

if the distribution of step heights, given by the probability 
density P(AE), is known. The process of linearizing will 
happen randomly throughout the array until, at some 
temperature T*, sufficiently many junctions have been 
replaced by Ohmic conductors that a continuous path 
involving only such conductors spans the array. At this 
point, the overall response will necessarily also be linear 
and the threshold must have reached zero: Vt(T*) = 0. 

An upper limit on when this point is reached can be 
obtained from percolation theory by considering the two 
classes of junctions as two types of bonds between neigh- 
boring sites. At small overall bias, we can label the non- 
linear junctions as insulators and the Ohmic ones as con- 
ductors. If a (temperature-dependent) fraction p{T) of 
all junctions in the array has been linearized, and in the 
absence of correlations between neighboring junctions, 
the first continuous path of linear conductors across the 
array occurs, on average, at a critical fraction pc- Here 
Pc is the bond percolation threshold which depends only 
on lattice topology and dimension (for corrections due to 
correlations see Section IV). Using Eq. 12 we thus find 
T* through 

piT*)=Pc. (3) 

As a consequence of these considerations, the global 
threshold will be a decreasing function of temperature 
and approach zero as p — > pc- Hence, to first order, 

VtiT) = Vt{0)a^PiT)/pc). (4) 

In order to proceed and find the linearized fraction of 
junctions, p(T), we need to know more about the actual 
distribution P{AE) of energy costs. It will be calculated 
in detail in Section HI. However, a few important aspects 
are already clear from Eq. 13 In particular, since Pc/2 is 
no larger than 1/4 for 2D lattices, we have to integrate 
over only a small portion of P{AE) in order to reach a 
significant suppression of the threshold. If P{AE) docs 
not change much over this range, we find 

p{T) « 26fcsTP(0) (5) 

and p{T) is proportional to temperature. The relevant 
energy scale, 1/P(0), can be thought of as an effective 
charging energy, while b depends only on the shape of 
the internal energy distribution of the metal particle and 



thus is independent of topology, dimensionality and the 
effects of coupling. 

We will see in Section HI that this is a reasonable ap- 
proximation not only for the case of zero capacitive cou- 
pling, but even more so when nearest neighbor coupling 
is included. Physically this is so because coupling flat- 
tens out the polarization-induced disorder in the energy 
landscape and small energy costs become more probable 
so that P{AE) decays slower for small AE. Combining 
Eqs. ^andOwe see that the normalized threshold decays 
linearly with temperature according to 




= 1 - i.8kBTP{0)/pc, (6) 



where we have used the result b = 2.4 obtained earlier. 
In analogy with the T — result Eq. ^ the right hand 
side of this equation represents a(T)/a, the temperature- 
dependent number of up-steps in the optimal path nor- 
malized by the number at zero temperature. 

Equation^is a central result of this paper. It predicts 
a linear depression of the global threshold with tempera- 
ture, with a prefactor 2bkBP{Q)/Pc that is universal and 
does not depend on the details of the threshold distribu- 
tion. 



III. ENERGY COST DISTRIBUTION 
INCLUDING NEAREST-NEIGHBOR COUPLING 

To calculate P{AE), we start from the electrostatic 
energy of a system of capacitors, 

E=\Y.{q^ + Q^)CT^HqJ+QJ), (7) 

where the qi are quenched, offset charges and the Qi 
are mobile charges (equal to an integer multiple of e = 
— 1.6 X 10~^^C or zero). The C,^^ are elements of the 
inverse capacitance tensor. Note that C^^, in the stan- 
dard definition of the capacitance tensor, does include 
contributions from coupling to nearest neighbors if such 
couphng is present. 

We define the energy difference before/after tunneling 
of a single electron from site 1 to site 2 as 

AE ^ Eq^^o^Q^^e - EQ^=e,Q2=o- (8) 

In the absence of any quenched charge disorder {qi — 0) 
we have AE = 0, and there is no cost associated with 
moving charges around inside the array. In other words, 
there is no Coulomb blockade of tunneling (even though 
Afij > 0) and the current-voltage characteristic will be 
linear. 

Now imagine a fiat, random distribution of quenched 
polarization charges in the range qi e [— e/2, -|-e/2]. As 
before, this range suffices because larger offsets will be 
compensated by mobile charges of magnitude e. In the 
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FIG. 1: Ten-sphere subsystem in a triangular lattice. The 
electron transfer occurs between sites 1 and 2. The other 
sites are the nearest neighbors. 



limit of negligible capacitive coupling between sites con- 
sidered for now, this leads to 

AS = e (gi - (72) Cfi^ 

To deal with nearest neighbor capacitive coupling, we 
focus here on the case of a close-packed, triangular lattice 
simply for the sake of having a concrete picture in mind 
and for direct comparison with experiments. In general, 
any lattice type can be treated the same way and the 
differences affect only the quantitative results for the ca- 
pacitance tensor elements. 

We consider a subset of the triangular lattice consisting 
of 10 spheres: two central sites (#1 and #2) participating 
in the tunneling event and their 8 surrounding neighbors 
as in Fig. ^ Keeping only nearest neighbor elements and 
taking Qj = for j > 2, 

AE = e (gi - (Cr/ - 0^2') + 

eCi2 (93 + 94 + 95 - qi - qs - 99) • 
Defining 7 = Cj^j^/Cjj^'^, we write AE as 



AE = e'Cil{[l-j]iq^ 



'72 j 



7 + 94 + 95 - 97 - 98 - 99)}- 



(9) 



The terms in round brackets, containing the qi, are sums 
of 2 or 6 random variables. The maximum value for AE 
is achieved if the appropriate limiting values (+e/2 or 
— e/2) are inserted for the qi. This gives 

AE,na. = e^C^,^ (1 + 27) . 

Without capacitive coupling to neighbors, AEmax can 
be written as AE^ax = e'^/Co, where Co = Aneeor, is the 
capacitance of a single sphere of radius r embedded in a 
medium of dielectric constant e. The key points emerging 



from equations |S1 and |51 are that the system energy cost 
associated with a tunnel event is not equivalent to the 
change in chemical potential of a single site, and that 
existence of a range of polarization charges qi gives rise 
to a distribution of energy costs AE. 

To calculate the full distribution P{AE) of energy dif- 
ferences, we need to first find the distributions P2{x) and 
Peix) resulting from the addition of 2 or 6 random vari- 
ables. In general, the probability of obtaining a value 
X = xi + X2 + ■■■ + Xn from the sum (or difference) of n 
independent random numbers Xi can be calculated from 
their recursion relation: 



Pn{x) 



dX'Pr,^l{x~x')Pi{x'). 



Using Fourier transform to convert the convolution 
into a product, we get Pn{£,) = Pn-i(C)^i(0- This leads 
to F„(0 = P{\£) = [sin(e/2)/(e/2)]", or 

2 1"°° sin" £ 
Pn{x)^- — ^cos(2fx)de 

TT Jo ? 

Specifically, for n — 2 and n — 6 this integral can be 
solved analytically and gives 

P^{x) = (|a; - 1| + |a; + 1| - 2|a;|)/2 
Pq{x) = (|a;-3|^ + |a; + 3|^-6|x-2|^-6|a; + 2|'V 
15|a; - 1|^ + 15|a; + 1|'^ - 20|x|^)/240 

The probability distribution of AE in Eq. |^ is then 
given by 



P{AE) 



1 



+ 00 



1 



-P2 



AE - AE' 

AE' 



dAE'. 



(10) 



The shape of this P{AE) is triangular with apex at 
AE = 0. Depending on 7, the shape is rounded near the 
top ( where AE' —>■ 0) and curved outward near the bot- 
tom (as AEmax is approached). The amount of round- 
ing/curving increases with 7 (Fig. Specifically, for 
negligible coupling (7 = 0), P{AE) becomes the distri- 
bution of differences between two random variables 



P{AE) = 1/AEr, 



\AE\/ {AEmax) 



(11) 



This is a simple triangle with P(0) ~ \/AE„iax and 
base extending from —AE,nax to + AEmax- Fig- [21 
shows P{e) as a function of the normalized energy cost, 
e^AE/{e^C^^). 

Using Eqs. OlandlTTlfor 7 = 0, we find that the fraction 
of linearized junctions is 



P{T) = 



2bkBT 



hkeT 



(12) 



For a 2D triangular lattice Pc — 0.347 so that the tem- 
perature at which an Ohmic conducting path percolates 
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of the capacitance tensor elements C^^ and • In 
Section IV we present numerical results for a range of 
coupling strengths and show how these tensor elements 
depend on the ratio of center-to-center distance, L, to 
particle radius, r. As particles get closer and L/r 2.4, 
7 reaches 0.4 and the approximation leading to Eq. 
breaks down (see also Fig. below). Furthermore, for 
very large interparticle coupling, next-nearest neighbor 
interactions will become significant and correlations be- 
tween energy-steps may become more important. 

We can repeat the above derivation of P(0) for a one- 
dimensional linear chain of particles. In this case, we 
consider 4 sites in a row with an electron moving between 
the two central sites. Now P{AE) contains the integral of 
a product of two P2 functions. We find that for 7 < 1/3, 



FIG. 2: Probability distribution of the energy cost of tun- 
neling between sites 1 and 2 in Fig. The distribution of 
£ = AE/{e'^C^-^) is plotted where AE is the change in the 
system energy due to tunneling. C^^^ is the diagonal element 



of the inverse capacitance matrix and 7 = Cj^j^/C 



1 

11 • 



across the lattice, defined in Eq. El by p{T*) — pc, is 
reached at bhsT* / AE^ax = 0.192. This value is small 
enough that, to good approximation, Eq. |S1 holds and 
the quadratic term in Eq. ^1 can be neglected for all 

For finite capacitive coupling between nearest neigh- 
bors, 7 > 0, P{e) in Eq. ^| can be expanded around 
e = to obtain 



P{e) = P(0) 



0.55 



7(1-7)^ 



(13) 



The linear term disappears because the distribution has 
a rounded top near e = (Fig. Consequently, correc- 
tions to Eq. |31are of order {bkBT* / AEj^axY and thus 
smaller than in the case of zero coupling (see Section IV 
for numerical integration results for T*). Therefore, the 
linear decrease of Vt with temperature in Eq. |21 holds to 
even better approximation. The first term, P(0), in Eq. 
1131 can be found straightforwardly from the integral in 
Eq. ^|as long as 7 is sufficiently small. This leads to 



P(0) = 



and finally. 



1 



1 2'y 
77 I xPp,{x)dx 

1-7 (1-7)^/0 ^ ' 



P(0) 



1 1 - 1.577 

e^cri^ (1 - 7)' ■ 



(14) 



Note that P(0) depends only on the geometry of the 
system and is independent of all details of the quenched 
charges (as long as they can be assumed uniformly ran- 
dom). This allows us to obtain P(0) from calculations 



P(0) 



ID 



I 1 - 47/3 



(15) 



One final aspect concerns how the zero-temperature 
threshold Vt (0) in Eq. is affected by capacitive cou- 
pling between neighboring particles. In MW's argument 
leading to Eq. ^for the uncoupled case, the factor e/Cg 
came from an increase in local potential corresponding to 
one full electronic charge. With capacitive coupling, the 
increase in local potential due to an electronic charge will 
be less as it effectively spreads out over the neighbors. 

In order to reach the threshold for conduction, we still 
have to add approximately one electron to the array for 
each up-step in a path. To first order, the average local 
change in potential associated with adding an electron 
is eCii, where decreases with increasing coupling. 
As before, a is the number of up-steps in the optimal 
path at T = divided by the length of the array. The 
optimal path is the one with the fewest number of up- 
steps. Let us define Vq as the average increase in external 
bias required to overcome an up-step. We then can think 
of the voltage threshold as a product of two quantities: 
the number of up-steps {aN) and the cost in bias per 
up-step (Vb ~ ^^ii)- Modifying Eq. ^ we are led to 



Vt{Q) = aNVo 



(16) 



Note, however, that this relation is only an approxima- 
tion and that a full calculation is a formidable problem 
for 7 > 0. The reason is that now local changes in poten- 
tial depend strongly on the quenched charge configura- 
tion as well as on other mobile charges arriving on nearby 
particles. In 2D, in particular, this complex interaction 
poses a challenge not only for analytical calculations but 
also for simulations. On the other hand, ID simulations 
can be carried out straightforwardly and can be used to 
gauge the validity of Eq. ^1 This will be done in the 
next section. 
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IV. NUMERICAL CALCULATIONS AND 
CHECKS 

In order to use Eqs. Eland^l we need to know certain 
elements of the inverse capacitance matrix as well as the 
value of a appropriate for a given lattice. Both of these 
can be obtained from numerical calculations as we detail 
in this Section. In addition, simulations allow us to per- 
form a number of checks of the assumptions underlying 
the model developed in Section III and they provide a di- 
rect test for the effect of correlations that were neglected 
in its derivation. In the following figures, we normal- 
ize capacitances by the capacitance, Co, of an isolated 
sphere and energies by e^/Cg, the maximum energy cost 
for tunneling between capacitively uncoupled particles. 

Inverse Capacitance Matrix. To calculate the capaci- 
tance matrix of the 10-sphere system in Fig. ^ we used 
FASTCAP, a capacitance extraction program developed 
at MITi2£ The program implements a preconditioned, 
adaptive, multipole-accelerated 3D capacitance extrac- 
tion algorithm developed by Nabors et a/..-^ Each site 
in the system was represented by a spherical, 1200-panel 
polygon. Center-to-center distances between 2.1 and 20 
times the radius were examined. (For L/r — 20, we used 
a 104-panel sphere approximation so as to not run out 
of computer memory.) The output of the program is a 
10x10 capacitance matrix C in units of pF for spheres of 
radius Im. We then inverted this matrix in Mathematica 
to find C~^. Since capacitance is directly proportional to 
the scale of the system, and to the dielectric constant, we 
can remove these dependences by scaling all capacitance 
elements by the self-capacitance of an isolated sphere. 
We will do this in all the figures to give a general result. 

Fig. O shows the effect of coupling on the 1-1 and 1-2 
elements of the inverse capacitance matrix. Note that 
the self-capacitance Cn and thus depends on in- 
terparticle coupling because nearby spheres can polarize 
when a charge is added to the central sphere, decreasing 
the overall energy cost of the charge addition. However, 
as Fig. 121 shows, for values oi L/r > 3 the change in 

due to nearest-neighbor coupling is small, and 
remains within 10% of I/Cq. Typical experimental val- 
ues for close-packed, dodecanethiol-coated 6nm particles 
give values L/r of about 2.74^ As Fig. |21shows, the off- 
diagonal element depends less strongly on L/r than 
the diagonal element. Thus, the increase in 7 with de- 
creasing L/r below a value of about 3 is largely due to 

'-'11 ■ 

We also note that the interparticle capacitance C12 
depends on having extra neighbors. For example, for 
L/r — 2.67, C12 in the 10-sphere system of Fig. ^is only 
71% of the value obtained for two isolated spheres. Thus, 
it is essential to look at the system as a whole and not to 
assume isolated spheres. In order to check whether or not 
the 10-sphere system is sufficient, we added another ring 
of spheres to Fig. 1, creating a 24-particle subset of the 
triangular array. We then calculated the full capacitance 
matrix for the 24-particle system. For L/r = 2.1, we 
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FIG. 3: EflFect of coupling on the elements of the inverse 
capacitance matrix for a 10-particle triangular system. Cou- 
pling increases as the center-to-center spacing, L, normalized 
by the radius, r, decreases. As coupling increases, the in- 
verse self-capacitance, Cii , decreases and the inverse inter- 
particle capacitance, C^j^, increases. We normalize by the 
self-capacitance of an isolated sphere, Co, to assure that the 
values plotted are independent of r and the dielectric con- 
stant. 



found that the changes in Cn and C12 are less than 
1%. Consequently, we take the 10-sphere system as a 
sufficiently good approximation to the triangular array. 

The results of the capacitance matrix calculation are 
in contrast to approximations^ which estimate the effect 
of capacitive coupling by adding to the capacitance of 
an isolated single particle, Co, the interparticle capaci- 
tance C12 for each neighbor. In particular, over the range 
2.1 < L/r < 4 the estimate Cn « Co + 6C12 for a trian- 
gular lattice gives about twice the value for Cn obtained 
numerically using FASTCAP. 

In principle, FASTCAP will give the full capacitance ma- 
trix of the 10-particle system in Fig. and thus take into 
account several longer range couplings. However, here we 
limit the discussion to nearest neighbor coupling. First, 
we examine the zero-temperature limit. 

Conduction Threshold at T = 0. To calculate numeri- 
cally the onset of conduction at T = 0, we follow MW's 
model and look for paths across the array that minimize 
the number of up-steps. For this, we use a variant of 
the well-known Dijkstra optimal path finding algorithm, 
the "bottleneck algorithm."— For each site, we define an 
offset charge qi. If qi > qj, then i-to-j is considered an 
energy up-step in the uncoupled case. While we cannot 
use this method to find the full current-voltage charac- 
teristics, it provides a very fast and effective way of de- 
termining the validity of Eq. ^and it allows us to extract 
the geometrical prefactor a. As defined in Section III, a 
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FIG. 4: Charge front in a 2D triangular lattice as a function 
of external bias. The mobile charges, shown in dark gray, are 
able to penetrate further into the array from a reservoir on 
the left as bias is increased from left to right in the 3 pictures. 
The simulations on a 100x100 array were carried out using 
the "bottleneck" algorithm from Ref. 22. 



is the number of up-steps in the optimal path at T = 
divided by the length of the array. Note that our defini- 
tion of a differs from MW who define a = Vt(0)Co/iVe. 
The two definitions only agree in the uncoupled case. 

We can also numerically obtain the charge front as it 
propagates across the array for voltages below threshold. 
To do so, we find all the sites that can be reached in less 
than a given number of energy up-steps. In Fig. ^ we 
show three snapshots from a simulation on a triangular 
lattice with increasing bias from left to right. The ad- 
vancing charge front is seen as the right-hand edge of the 
dark gray region. 

Let us first consider the uncoupled case. For all types 
of lattices investigated, we find that Vt(0) increases lin- 
early with N ^ as predicted by Eq. ^] For a 2D 
square lattice MW reported a = 0.338(1) using Monte 
Carlo simulations. The bottleneck algorithm gives a — 
0.329(7) for a 160x160 square lattice (averaged over 1000 
trials). For honeycomb and triangular arrays (100x100 
array, 1000 trials) we find a = 0.301(9) and a = 0.226(8) 
respectively. 

What is the effect of coupling on the number of up- 
steps in the optimal path? A "step" A£' between two 
sites is not just (g^ — qj)/Ca, but now takes into account 
all neighbors, as in Eq. |^ However, since a does not de- 
pend on the magnitude of the up-steps, we do not expect 
a large effect. This is borne out by the simulations. In 
ID, a is not affected by coupling even for C12/C11 >> 1. 
In a 2D triangular array, we find that a depends only 
weakly on coupling. For L/r = 2.1, a decreases by about 
10% from its uncoupled value; for L/r >b, a has essen- 
tially the uncoupled value of 0.226. 

In order to compare our model more directly with lit- 
erature results for the global threshold in the coupled 
case, which are available only in IDji we simulated a 
ID chain of sites. In this simulation, we only con- 
sider self-capacitance iC^^) and nearest-neighbour ca- 
pacitance (Cj~2^). An electron moves forward from site i 
to i -I- 1 if Ai?i^i+i < 0, where IS.E is calculated from Eq. 
[3 considering both offset charges q and integral charges 
Q from all previous tunneling events on all sites. 

The external bias is raised in increments much smaller 
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FIG. 5: Average external voltage bias per up-step, Vb, at 
threshold in a ID chain of spheres as a function of interpar- 
ticle coupling at T = 0. The vertical axis is normalized by 
the bias per up-step in the uncoupled case, e/Co, where Co 
is the self-capacitance of an isolated sphere. The horizontal 
axis is the interparticle capacitance, C12, normalized by Co. 
The data from our ID simulation (open stars) are compared 
with simulation results from Ref. 1 (full squares) and two 
analytical approximations (open triangles and filled circles). 



than eC^i to inject electrons into the system. Electrons 
are allowed to propagate forward and rearrange to find 
the minimum energy state of the system before increasing 
the bias again. Vt is the external bias value for which 
the first electron reaches the far end of the chain. For 
each disorder realization in a 100-site chain, we count 
the number of up-steps and then raise the bias to find 
the threshold. As mentioned in the previous section, for 
finite coupling, there is no unique cost in bias per up- 
stcp, but rather a distribution. Fitting the average cost, 
Vo, to a quadratic function for 7 < 0.4, we find 



1 



1- 1.937+ 1.5372 + 0(7^ 



(17) 



Results from this simulation are shown in Fig. (SJ where 
we plot the average cost per up-step, Vb, normalized to 
the uncoupled value, as a function of C12/C0 in a ID 
chain. This is compared to the approximations Vb « 
eC^l from Eq. [Hand Vb « l/eP(0) using the ID result, 
Eq. ^|for P(0). Also shown are three data points from 
MW's Fig. based on a full simulation of the current- 
voltage characteristics of a chain. (Note that MW use 
a different normalization in their Fig. 1, i.e., they plot 
VtCo/eN, and extend the simulations to larger coupling 
strengths.) 

Conduction Threshold for T > 0. As a next step, we 
add temperature to the simulations. In the 2D algorithm 
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FIG. 6: Effect of temperature on an 2D triangular array with 
quenched charge disorder. As temperature is increased from 
left to right in the 3 pictures, mobile charges (in dark gray) 
can penetrate deeper into the array without energetic cost. 
When a percolating dark grey path spans the array from left 
to right, the global threshold bias for conduction reaches zero. 
The simulations on a 100x100 triangular array were carried 
out using the "bottleneck" algorithm from Ref. 22. 
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that finds the optimal path across the array we have di- 
rect access to all bonds, and thus energy costs, AE, for 
moving an electron between any pair of neighboring sites. 
This allows us to test the validity of the linear decrease 
of the threshold with increasing temperature predicted 
by Eq. |B1 As temperature increases, steps with mag- 
nitudes smaller than a threshold energy, AEth — bksT 
are thermally erased. In the path-finding algorithm, we 
count all steps less than AEth as "down-steps", that is, 
they do not cost any energy. We then traverse the energy 
landscape to find the least-cost path as before. Fig. El 
shows three snapshots from the simulation. The thresh- 
old AEth, and thus temperature, is increased from left to 
right. In dark grey we show all sites reachable without 
cost from the left edge. 

In Fig. [3 we plot a{T) as a function of the effective 
temperature, AEth, for various degrees of coupling. In all 
cases, we find that the number of up-steps in the optimal 
path decreases with increasing threshold energy approx- 
imately linearly. (The deviations from a strictly linear 
decrease, close to a{T) = 0, come from a finite size effect: 
the simulated arrays contained 100x100 sites, so around 
a(T) = 0.01, the average number of up-steps reaches 1, 
below which the average is fractional and thus no longer a 
physical measure.) We also see that for a given "temper- 
ature" the threshold decreases with increasing coupling. 
Furthermore, the temperature T* at which a{T) = de- 
creases with increasing coupling (see also Fig. ^jfl). In 
accordance with the results in Fig. [SJ these trends are 
most pronounced for small L/r and saturate near the 
uncoupled behavior for L/r > 5. 

Percolation and Correlations. In the analytic calcu- 
lation of T* in Section II, we found the fraction of lin- 
earized bonds, p{T), through Eq. [3 and defined T* as 
p(T*) — pc, where Pc is the bond percolation threshold 
in the lattice under consideration. This procedure relies 
on two assumptions that we now test. 

First, the basic idea of our tunneling model is that none 
of the down-steps cost energy. Implicit in Eq. |21is a some- 
what more restrictive criterion, namely |A£'| < fofc^T, 



FIG. 7: Decrease of the number of energy up-steps in 
the optimal path with temperature. ce{T) is defined as the 
temperature-dependent number of up-steps in the least cost 
path across the array divided by the array length. The ef- 
fect of thermal fluctuations is introduced by counting up- 
steps only if they exceeded a cut-off energy AEth = hknT. 
Data are shown from simulations on a 2D triangular lattice 
containing 100x100 spherical particles for three different cou- 
pling strengths, parameterized by the ratio of center-to-center 
distance, L, to sphere radius, r. The inset shows the col- 
lapse of the curves upon normalization by a = a(T = 0), 
and by 1/P(0), the relevant energy scale for temperature- 
dependence. 

requiring that the path be along only those bonds (cor- 
responding to either up- or down-steps) that had been 
linearized by thermal fluctuations. There may be en- 
ergetically much more optimal, but more asymmetric, 
paths that take advantage of those larger-energy down- 
steps that have not yet been linearized. In this situation, 
we are starting with a lattice with all down-steps in place 
and ask when the system-spanning path forms in the pro- 
cess of adding up-steps of increasing size. 

Second, setting p{T*) = pc and using literature values 
for Pc assumes that the usual rules for bond percolation 
apply and, in particular, all bonds are placed completely 
randomly. However, while the site energies are from a flat 
distribution, the energy differences between sites are cor- 
related. For example, in the triangular lattice in Fig. ^ 
the energy differences between sites 1 and 2 and between 
sites 2 and 6 completely specify AE between sites 1 and 
6. Therefore, Pc may not necessarily provide an accurate 
value for the threshold. Finally, even in the absence of 
correlations, the finite size of arrays corresponding to ex- 
perimental situations (with N no more than a few 100) 
may lead to a small correction to pc as listed for infinite 
lattices. 

We investigated these questions for a variety of lattice 
types by calculating numerically the threshold pa , defined 
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0.970(7) 


0.968(39) 


0.945(45) 
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0.653 


0.643(9) 


0.627(15) 


0.539(15) 


4 


0.5 
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0.384(10) 
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0.292(6) 
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0.250(5) 


0.307(7) 
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TABLE I: Percolation coefficients for different coordination 
numbers z, calculated for 200x200 arrays and averaged over 
200 trials, ps is the average fraction of bonds that need to 
be linearized in the whole array such that the first system- 
spanning path appears containing only linearized bonds (both 
up- and down-steps). Pa is the average fraction of bonds lin- 
earized in the array for the first system-spanning path con- 
taining non-linearized down-steps as long as all up-steps are 
linearized, ps and pa are for systems with random site ener- 
gies. Pc is the bond percolation fraction for the uncorrelated 
bond percolation. The theoretical values Pc,th are taken from 
Ref. 19 and presented for comparison to indicate the extent 
of finite-size effects. 



as the average fraction of bonds required for percolation 
under the asymmetric condition AE < bksT*, and the 
threshold ps, defined as the average fraction of bonds 
required for percolation under the symmetric condition 
|A£;| < bkBT*. Both are hsted in Table 1, together 
with Pc as obtained from the same lattice but with ran- 
domly assigned bond energies rather than random site 
energies, z is the coordination number, the number of 
nearest neighbors of each site. The lattice with z = 2 
consists of 200 parallel ID wires. 

In Fig. |S1 a comparison between ps and Pc gives a sense 
of the relevance of correlations which increase Ps roughly 
linearly with increasing z (ignoring the case of z = 2). 
At the same time, antisymmetric paths involving large 
down-steps give a threshold fraction, pa, that is system- 
atically lower than ps by about 15%. Intriguingly, and 
quite unexpectedly, for lattices with z = 6 to z = 8 the 
contributions from correlations and asymmetry appear 
to cancel each other to a large extent so that pc provides 
an excellent estimate of the "true" value, Pa- Thus, using 
Pc in Eq. Oto estimate T* should give very reasonable 
estimates for experiments on self-assembled nanoparticle 
layers. The small difference between the theoretical pc 
and the value from the simulation shows the insignifi- 
cance of finite size effects for 200x200 arrays. 

Distribution of Energy Costs. The last approximation 
in the model we wish to test is the replacement of the 
integral in Eq. Hwith pc = iMsT* P{0). To find the 
full distribution of energy costs for the nearest neighbor- 
coupled 10-sphere system shown in Fig. ^ we used a 
Monte Carlo routine. Offset charges from a uniform ran- 
dom distribution [— e/2, -|-e/2] were assigned to each of 
the 10 sites. Using the capacitance matrix as calculated 
from FASTCAP, the energy cost AE associated with tun- 
neling from site 1 to site 2 was found from Eq. |51 for 
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FIG. 8: Comparison of symmetric and asymmetric percola- 
tion conditions for various 2D arrays with coordination num- 
bers, z. Ps, Pa and Pc are defined in Table 1. All simulation 
data are for arrays of size 200x200 and averaged over 200 
disorder realizations. 



each disorder realization. P{AE) was then obtained from 
sampling AE for 10^ offset charge realizations for each 
value of L/r. 

Fig. 1^ shows the normalized peak probability density 
P(0)e^/Co as a function oi L/r. P(0)e^/Co only depends 
on L/r since both Co and 1/P(0) are proportional to r 
and e. We compare the simulation value with the ap- 
proximation in Eq. ^| Knowing the full distribution 
from Monte Carlo simulations allows us to find, without 
approximations, the critical temperature T*, where the 
voltage threshold goes to zero. According to Eqs. |21and|21 
this is done by integrating P{AE) out to the point where 
the area under the graph corresponds to pc. In Fig. |Ht>, 
we compare the results of numerical integration with the 
analytical approximation T* = pc/(2&fcsP(0)) (Eq. I^J 
with P(0) from Eq. E]for a 2D triangular lattice. 



CURRENT- VOLTAGE CHARACTERISTICS 
ABOVE r* 



Next, we investigate the behavior for temperatures 
T ~ T* and above. Within our model, T* is defined 
as the temperature at which there are just enough local 
junctions linearized to span the array at zero bias and re- 
move the global threshold. In other words, with increas- 
ing temperature the nonlinear current-voltage (/ — V) 
characteristics, described by the powerlaw / ~ (V" — 
Vf(T))^, have been linearly shifted to the left until, at 
T*, they first reach the origin with finite slope. This gives 
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FIG. 9: Effect of capacitive coupling on the energy scales 
associated with the temperature-dependence of the global 
threshold for conduction. L is the center-to-center distance 
between neighboring spherical particles and r is their radius, 
a) Changes in the peak of the energy-cost distribution, P{0), 
as a function of coupling. 1/P(0) plays the role of an ef- 
fective charging energy; in this plot it has been normalized 
by the maximum energy cost in the uncoupled case, e^/Co- 
Results from Monte Carlo data (full squares) and from the 
approximation given by Eq. 1141 (open circles) are shown, b) 
Changes in T* , the temperature at which the voltage thresh- 
old of the array becomes zero, for a triangular lattice (perco- 
lation threshold, pc = 0.347). T* has been normalized by Tq , 
the value in the uncoupled case. Closed square symbols are 
data from numerical integration of the energy cost distribu- 
tion, P{AE), as obtained by Monte Carlo simulation; open 
circles show the approximation given by Eq. |^ 



rise to a finite zero-bias conductance, go = dl /dV\v=o- 

For T > T*, additional linearized junctions provide 
parallel paths across the array and the zero-bias conduc- 
tance increases. For any given temperature, however, 
an increase in bias will eventually provide sufficiently 
high local voltage drops to involve portions of the array 
with junctions not yet linearized by thermal fluctuations. 
Thus, at sufficiently high bias, the I — V characteristics 
will change back from Ohmic to the original nonlinear 
powerlaw behavior with temperature-independent expo- 
nent C- 

These considerations correspond to a picture in which 
the nonlinear I — V curve of a 2D array emerges from 
summing the contributions from all paths carrying cur- 
rent at a given bias voltage (this is the essence of the 
scaling derived by MW). As the bias is increased above 
Vf(T), more and more paths will open up as their up- 
steps are overcome. Above T*, when the global threshold 
disappears along the optimal path, there still are many 
other paths that have finite thresholds and are accessi- 
ble at higher bias. These thresholds will keep decreasing 
linearly with temperature as more and more steps are lin- 
earized by thermal fluctuations. As a consequence, the 
high-current, powerlaw portion of the I ~ V curves will 
shift to lower and lower bias voltages. 

We therefore expect the I — V curves to collapse, by 



a simple horizontal shift, onto a master curve not only 
for T < T*, but, at least with their high-bias portion, 
also for T > T*. In the latter regime a true threshold for 
global conduction obviously no longer exists and Vt(T), 
as obtained from the shift required for a high-bias I — V 
collapse onto a common power law, should be thought of 
as an effective threshold value. This effective Vt{T) will 
be negative for T > T*. 

These predictions are borne out by experiments. As 
first observed by Ancona et al^ and also seen in Fig. ^1 
the linear decrease of Vt with T continues past T* and 
into a regime of negative effective threshold values. (Our 
simulation results in Fig. |3 are not based on calculating 
full I—V curves. Since our method finds the first system- 
spanning path, it is no longer applicable above T* where 
such paths exist already at zero applied bias.) 

We note that our model provides a natural explana- 
tion for this smooth cross-over from positive to negative 
Vt{T) which does not invoke additional mechanisms. In 
particular, it does not bring into play activation over the 
tunnel barriers for T > T*, as proposed by Ancona et 
al..^ Because such barriers are set by the properties of the 
alkanethiol ligands separating adjacent nanoparticles, we 
expect barrier heights in excess of several eV, correspond- 
ing to the first available electronic states in the ligands. 
This is significantly higher than either fc^T or the volt- 
age drop per up-step, making hopping over the barrier 
highly unlikely. 

How does go depend on temperature? To address this 
point, we revisit a simplifying approximation made in 
the model, namely the coarse-graining in which non- 
linearized junctions were assumed to be unaffected by 
thermal fluctuations and to exhibit zero conductance be- 
low threshold. In principle, of course, flnitc T will always 
induce some zero-bias conductance. For a single junction 
this zero-bias conductance exhibits activated behavior, 
i.e., is proportional to exp(— C//A:bT), where the activa- 
tion energy, U = AE, is the energy cost required to move 
a charge across the junctionii^ 

Sufficiently far below T*, there always will be several 
junctions along the optimal path with AE » ksT. 
Consequently, the overall zero-bias conductance will be 
exponentially suppressed to a level where 170 is well ap- 
proximated by the coarse-graining approximation. 

Once T* is approached, linearized junctions for the first 
time form a system-spanning path, go will be dominated 
by the relatively few bottleneck junctions with the largest 
activation energy in the path, AE ~ ksT* . The majority 
of paths around the bottlenecks would involve junctions 
with much larger AE which therefore could shunt the 
bottlenecks only insignificantlyi2i24 

Therefore, near T* the overall, zero-bias array con- 
ductance, go, will display activated behavior similar to a 
single junction with U = bksT* = pc/(2P(0)). As be- 
fore, the key point here is that the activation energy is 
not simply the charging energy of an isolated grain, but 
is connected to the optimal path across an energy land- 
scape established by the quenched charge disorder. For 
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2D triangular arrays Pc/2 ~ 0.17 so that U is approxi- 
mately 1/5 of the effective charging energy 1/P(0). 

Above T*, additional paths in which all up-steps have 
been thermally erased will span the array. These paral- 
lel paths will contribute to go and modify the behavior. 
Taking the number density, D{T), of such paths to be 
proportional to the percolation conductivity of the lin- 
earized subset of junctions above Pc, we have D{T) ^ 
{p{T) -pcY, where p{T) is given by Eq. Hand t w 1.3 
in 2D«iS The overall zero-bias conductance at tempera- 
ture T then follows from integration of D{T')go{T') be- 
tween T* and T . This will give a powerlaw correction to 
the simple activated behavior, but we expect the optimal 
path established at T* to continue to dominate since its 
bottlenecks have the lowest activation energy of the set. 

The experimental data in Ref . il2| for the zero-bias con- 
ductance in 2D arrays above T* is compatible with simple 
activated behavior, with values for U that are within a 
factor two oi pd {2 P {{))). However, since T* can reach 
lOOK or more in these arrays, the remaining interval up 
to room temperature simply is not large enough to pro- 
vide a stringent test of the model predictions. The simple 
activated form for was also observed by Black et al. 
in arrays of Co nanoparticlesi^ 



VI. DISCUSSION 

The model described in the preceding sections provides 
a physical picture for the role of thermal fluctuations as 
they affect the nonlinear transport properties in systems 
with random local thresholds for conduction. One key 
aspect is that such systems cannot be described by as- 
sociating a single, fixed energy cost with charge motion 
from site to site. Instead, it is important to consider the 
fact that energy costs depend on both the Fermi levels 
of the particles involved in the tunneling process, and on 
the surrounding charge environment. This is the essence 
of Eq. |S1 which gives the change in electrostatic energy 
of the system as a whole, and, in the nearest neighbor 
approximation considered here, leads to the probability 
density of energy costs shown in Fig. |2] 

In particular, the energetics of the system are not 
determined simply by the change in the Fermi level 
(by some fixed single electron "charging energy" ) of the 
nanoparticle tunneled into. Energy costs derive from dif- 
ferences in the total system energy before and after a tun- 
neling event and thus involve at a minimum two particles 
or, with capacitive nearest neighbor coupling, ten parti- 
cles in a 2D triangular lattice (Fig. ^ or four particles 
along a ID line. In the absence of quenched charge dis- 
order, these energy costs vanish and mobile charges can 
move freely inside the array (the only costs are incurred 
when charges enter the edges of the array from one of the 
electrodes) . Therefore, the scaling of the global threshold 
with system size N in Eq. (the original MW result) or 
its modification for finite coupling strength, Eq. 1161 are 
a direct consequence of quenched charge disorder. 



However, as far as the global threshold is concerned, 
the detailed shape of the full energy cost distribution 
turns out to be not critical. Rather, as we showed by 
mapping the system onto an equivalent percolation prob- 
lem, the behavior is dominated by the lowest-cost per- 
centiles (up to Pc/2). To very good approximation, this 
is captured by the zero-cost, peak value of the probability 
density, P(0) (Eqs. El and El for 2D and ID systems, 
respectively). Increased capacitive coupling is found to 
have two effects: it flattens the energy landscape, thereby 
narrowing the width of the energy cost distribution and 
increasing the value of P(0), and it rounds off the peak 
of the distribution, making P{0) an even better approx- 
imation of the relevant portion of P{AE). 

There are several predictions that emerge from the an- 
alytic model. 

- First, the model predicts a linear decrease of the over- 
all, global threshold with temperature (Eq. EJ. This de- 
crease is directly proportional to fcBTP(O), the strength 
of thermal fluctuations measured relative to the effective 
charging energy 1/P(0). What is particularly appealing 
about this result is that all details about capacitive cou- 
pling and particle geometry enter through P(0), while 
the numerical prefactor, 4.8/pc, captures the underlying 
network topology and dimension through the percolation 
threshold, pc. Therefore, data from systems with similar 
structure but different particle sizes, spacings, or dielec- 
tric constants should collapse onto a "universal" curve 
when plotted as Vt{T)/Vt{0) versus ksTPiO). 

- Second, for sufficiently broad distributions of 
quenched charge disorder we expect that thermal fluc- 
tuations do not alter the basic character and roughness 
of the energy landscape. Thus we expect MW's zero- 
temperature powerlaw scaling of the nonlinear current- 
voltage (/ — V) characteristics, I ^ {V — Vf (0))^, to sur- 
vive at finite T with fixed exponent but with Vf(0) re- 
placed by Vt{T). In other words, the shape of the I — V 
characteristics remains unaffected by temperature while 
they are shifted linearly towards smaller threshold values. 

- Third, the threshold is expected to vanish and the 
nonlinear, Coulomb-Blockade-type current voltage char- 
acteristics are expected to change to linear, Ohmic be- 
havior near zero bias once the temperature exceeds T* = 
Pc/{4:.8kBP{0)) (Eqs. Eland 13). 

- Fourth, also for finite capacitive coupling, the zero- 
temperature global threshold value, Vt(0), can be written 
as a product of two quantities: the average number of 
up-steps encountered along the optimal path across the 
array, aN, which depends mainly on array geometry, and 
the average applied voltage per up-step, Vq, at threshold. 
We argued that, at least to first order, capacitive coupling 
can be taken into account by Vq ~ ^C^-^ (Eq. Il()|) instead 
of e/Co for the uncoupled case (Eq. 0). 

The robustness of the linear decrease of the global 
threshold for conduction with temperature is underscored 
by recent simulations of the full current-voltage charac- 
teristics for mobile charges hopping between traps in a 2D 
lattice. The charges interact via a long-range Coulomb 
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term and the trap depth at each site is chosen from 
a Gaussian distribution. In this work Reichhardt and 
Olson Reichhardt not only find behavior qualitatively 
similar to our model, but also track as a function of 
temperature the charge flow patterns beyond threshold. 
Their results support a basic, underlying tenet of our ap- 
proach, namely that small thermal fluctuations simply 
shift Vt(T) to smaller values while preserving the rough- 
ness of the energy landscape and thus the shape of the 
current- voltage characteristics. 

Conversely, such linear shift might be taken as an 
indicator of a sufficiently wide distribution of trapping 
depths or local threshold values: simulations by Rendell 
et al. of small 2D tunnel junction networks using the 
full, temperature-dependent "orthodox" Coulomb block- 
ade model show this shape-preserving shift cnerging once 
quenched charge disorder is introducediS 

The fact that P{0)e'^/Co, through C^^^ and C^^, 
depends only on L/r makes comparison with experi- 
ments straightforward as long as the array geometry 
is known. In particular, if L and r can be deter- 
mined from transmission electron micrographs. Fig. |2| 
together with Co — Aneenr give direct access to P{0) 
for triangular arrays, and thus make it possible to plot 
the normalized threshold, Vt{T) /Vt{0) as a function of 
scaled temperature, kBTP{0). Fig. [TUI shows such plot 
for temperature-dependent threshold data obtained from 
self-assembled, close-packed gold nanocrystal monolayers 
covering a range of array lengths (27 < N < 170) and ef- 
fective charging energies 96meV < 1/P(0) < 302meV.^'^ 
The threshold values were obtained from linear shifts of 
the full I — V curves onto a single powerlaw mastercurve 
for each sample (with temperature-independent exponent 
C = 2.25 ±0.1). 

All of these data are seen to cluster around the lin- 
ear decay with slope — 4.8/pc = —13.8 and x-intercept 
kBT*P{0) = pc/4.8 = 0.07 predicted by the model. Note 
that this data collapse contains no free parameters once 
Vt{T) has been measured (over a wide enough range to 
extrapolate reliably to Vt(0)) and P{0) been obtained 
from the array geometry. If direct access to L and r is 
not possible, but TV can be estimated from the electrode 
spacing, P{0) can be estimated using Eq. ^| together 
with the appropriate value for a listed in Section IV. 

The linear suppression of Vt (T) with temperature was 
also observed by Ancona et al^ in experiments on 2D 
arrays of gold nanoparticles, and by Bezryadin, Wester- 
velt and Tinkham.^** in studies of ID carbon nanoparti- 
cle chains. While micrographs allowing for a determi- 
nation of L and r were not available from either exper- 
iment, Bezryadin et al. found that the voltage thresh- 
old decreased as Vt{T) « Vt{0) - NksT/e. The au- 
thors also give the radius of the carbon particles as 
well as N and Vt{0). From this information we esti- 
mate a cost per up-step of about 0.2e/Co. From Fig. 
|S1 and capacitance calculations of a 4-particle ID chain, 
we find P(0) = 1.66(Co/e2). With pc = 1 for a ID 
chain our model gives 14 (T) = VtiO) - 4.8Vt{0)P{0)kBT 
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FIG. 10: Decrease of the normalized voltage threshold as a 
function of effective temperature variable fcsrP(O) for 6nm 
diameter, close-packed gold nanoparticle monolayers. (Repro- 
duced from Ref. 12 with permission.) 



from Eq. The second term can be written as 
-4.8[0.5iV(0.2e/Co)](1.66Co/e2)fcBT, where the term in 
square brackets is Vt(0) and am = O.Si Using the exper- 
imental parameters given by Bezryadin et al., our theory 
predicts Vt{T) = Vt(0) - O.SNksT/e, which is close to 
the experimental temperature-dependence. 



VII. CONCLUSIONS 

Our model describes the effect of temperature on the 
global threshold for conduction through arrays with dis- 
tributed local thresholds due to quenched charge disor- 
der. It resolves two long-standing issues, namely how to 
extend the original, T ~ scaling approach by MW to 
T > and how to include capacitive coupling between 
nearest neighbors. 

One key finding is a robust, linear decrease of the 
global threshold with temperature, in excellent agree- 
ment with recent experiments on close-packed nanocrys- 
tal arrays. This explains the experimental finding that 
powerlaw I — V characteristics resulting from Coulomb 
blockade effects keep their nonlinear shape to remarkably 
high temperatures while simply being parallel-shifted as 
T is increased. The model further predicts the existence 
of a cross-over temperature T* , above which the low- 
voltage portion of the I — V characteristics changes and 
acquires a significant zero-bias conductance that exhibits 
simple activated behavior. 

A second key finding is the identification of 1/P(0) 
as the relevant, effective charging energy, extending ear- 
lier results that did not treat capacitive coupling. Our 
approach explicitly takes into account the fact that 
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quenched charge disorder leads to a distribution of en- 
ergy costs for tunnehng, even for otherwise perfect lat- 
tices of identical junctions. Finally, we present numerical 
calculations that allow one to extract the relevant capac- 
itance matrix elements as well as P{Q) from knowledge 
of interparticle spacing and particle diameter. 
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